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ABSTRACT 

We show that short-term perturbations among massive planets in multiple 
planet systems can result in radial velocity variations of the central star which 
differ substantially from velocity variations derived assuming the planets are exe- 
cuting independent Keplerian motions. We discuss two alternate fitting methods 
which can lead to an improved dynamical description of multiple planet sys- 
tems. In the first method, the osculating orbital elements are determined via a 
Levenberg-Marquardt minimization scheme driving an N-body integrator. The 
second method is an improved analytic model in which orbital elements such 
as the periods and longitudes of periastron are allowed to vary according to a 
simple model for resonant interactions between the planets. Both of these meth- 
ods can potentially determine the true masses for the planets by eliminating the 
sin i degeneracy inherent in fits that assume independent Keplerian motions. As 
more radial velocity data is accumulated from stars such as GJ 876, these meth- 
ods should allow for unambiguous determination of the planetary masses and 
relative inclinations. 

Subject headings: stars: planetary systems 


1. Introduction 

Several thousand nearby stars are now being surveyed for periodic radial velocity 
variations which indicate the presence of extrasolar planets (Marcy, Cochran & Mayor 
2000). In the past year, the pace of discovery has increased, and there are now nearly sixty 
extrasolar planets known. 

Recently, systems with more than one planet have been found, and four ( v Andromedae, 
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GJ 876, HD 83443, and HD 168443) are now known. GJ 876 (Marcy et al 2001) provides 
an especially interesting case. In this system, a combined, two-Keplerian fit to the radial 
velocity data (see Table 1), suggests that the star is accompanied by two planets on orbits 
having a nearly commensurate 2:1 period ratio. The amplitudes of the star s radial velocity 
variations suggest minimum masses of 0.56 Mjup for the inner planet, and 1.89 Mjup for 
the outer planet. GJ 876, an M dwarf star with an estimated mass of 0.32 ±O.O5M 0 (Marcy 
et al 2001), is the lowest mass star known to harbor planets. 

For these orbital parameters, the mutual perturbations of the two planets in the 
system are considerable. Using a Bulirsch-Stoer integrator (with a timestep accuracy of 
AE/E = l.Oe — 15) we have computed radial velocity curves using the Table 1 orbital 
elements given in Marcy et al (2000) as an initial condition. For co-planar Keplerian orbits, 
the orbital elements Pi t 2 , ei t 2 , T\ , 2 , ^ 2 , an d are constant, and these, along with the 
mass of the star, serve to completely determine the positions and velocities of the planets 
at any time. When the interactions between the planets are non-negligible, however, the 
orbital elements change continuously, and so one must also specify an initial epoch in 
order to determine the motion at future times. Every starting epoch corresponds to a 
different resultant three-body motion. In figure 1, we demonstrate this effect by mapping 
the elements reported by Marcy et al. (2001) onto initial conditions corresponding to JD 
2450106.2 (the reported time of perihelion passage for the outer planet). The red line 
shows the radial velocity curve of the star which results from the superposition of the two 
Keplerian reflex motions. The black line shows the radial velocity curve resulting from the 
full three-body integration. After three orbits of the outer planet, the motion begins to 
deviate noticeably from the dual-Keplerian approximation. After several years, the motions 
have diverged completely. 

The rest of this paper is organized as follows: in §2 we show that self-consistent radial 
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velocity curves are required for systems such as GJ 876. In §3, we derive improved, fully 
self-consistent dynamical fits to the observed radial velocities of the GJ 876 system. In §4, 
we show how dual-Keplerian fits can be improved using an approximate analytic model for 
the interactions between two massive planets in resonant systems. Further applications are 
discussed in §5, which also serves as a conclusion. 


2. Dual Keplerian vs. Self-Consistent Fits 

The orbital elements given in Table 1 (taken from Marcy et al. 2001) were derived 
under the assumption that they are constants of the motion. However, for a system such 
as GJ 876, where the mutual planetary interactions are strong, the elements will change 
quite rapidly on observable timescales. We can therefore regard the parameters in Table 1 
as a set of osculating elements, which, given a particular starting epoch, correspond to a 
uniquely determined initial condition. 

Even with the assumption that sini =1 for both planets, the variety of motion 
corresponding to the starting conditions given in Table 1 is very broad. For some starting 
epochs, the planets are not in the 2:1 resonance, and the system experiences severe 
dynamical instabilities within five years. For other starting epochs, the planets undergo 
librations about the resonance, and the system is stable over timescales of at least 70 million 
years (Marcy et al 2001). One can thus ask the question: are there any starting epochs for 
which the osculating elements in Table 1 generate an evolution which is consistent with the 
observed reflex velocity of the star? 

We have computed synthetic radial velocity curves resulting from these osculating 
elements using 10,000 different initial starting epochs spaced one day apart. In each case, a 
uniform velocity offset was applied to the synthetic radial velocity curve in order to match 
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the first radial velocity point obtained by Marcy et al (2001) with the Keck telescope (t — JD 
2450602.1, v = 343.72ms -1 ). We then compared each synthetic curve to the remaining 
53 radial velocity observations obtained at Keck, and computed a reduced x 2 statistic for 
the fit. The best fit occurred for an integration corresponding to a starting epoch of JD 
2450671.98. The reduced x 2 value for this fit is 17.27, and the rms scatter of the velocities 
about the curve is 83 ms -1 . The best fit curve, along with the observed data, is shown 
in Figure 2. Given that the observational errors lie in the range 3-5 ms -1 , this degree of 
scatter indicates that the best dual-Keplerian fit is a poor match to the observed velocities 
when mutual planetary perturbations are taken into account. 


3. A Self-consistent N-body Minimization Scheme 

The experiment described above demonstrates that it is essential to include mutual 
planetary perturbations when making fits to velocity observations of planetary systems 
resembling GJ 876. One way to do this is to attempt a fully self-consistent fit which employs 
N-body integrations to produce a synthetic reflex velocity curve for the central star. Starting 
with the best dual-Keplerian fit to the Keck data, we have used a Levenberg-Marquardt 
algorithm (Press et al 1992) to iterate an improvement to the osculating orbital elements 
reported in Table 1. Our implementation of the algorithm examines how the x 2 value of the 
fit depends on variations of all 10 orbital elements, and attempts to find a set of elements 
for which the reduced x 2 fit is at a minimum. The N-body integrations were done using the 
Bulirsch-Stoer integration package developed by Laughlin Sc Adams (1999). A preliminary 
investigation of systems in which sin i = 1 for both planets has shown that there are many 
isolated minima within the ten-dimensional parameter space associated with every starting 
epoch. 

Using the Levenberg-Marquardt algorithm, we have found a self-consistent model for 
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the radial velocity data which has a reduced x 2 value of 1.93 and an rms scatter of 12.0 
ms -1 . This fit is shown in Figure 3, and represents a large improvement over the fit shown 
in Figure 2. The osculating elements are given in Table 2. We would like to stress that this 
fit was the result of a preliminary investigation, and that better fits can almost certainly be 
found (even for sini=l). The fact that the best dual-Keplerian fit of Table 1 has a similarly 
low reduced x 2 value, despite diverging substantially from true three-body motion, indicates 
that the present data set can be well-modeled by a variety of functions. It is likely that more 
observations will be required in order to secure all 14 orbital elements (allowing for mutually 
inclined planetary orbits). It seems clear, however/that as more data are acquired, the true 
masses of the planets will be revealed. The mutual perturbations between the planets will 
provide additional information which overcomes the degeneracy which previously made it 
impossible to determine the inclinations (and thus the true masses) of the planets. 


4. Improved Analytic Approximations 

The foregoing Levenberg-Marquardt minimization method is a potentially powerful 
technique for determinining all of the orbital parameters of the system in a completely 
self-consistent fashion. However, in the absence of a good initial model for iteration, 
it is difficult to locate the global minimum for the system. As we have shown above, 
the best dual-Keplerian fit provides a very reasonable starting model for the Keck data 
alone. However, the best dual-Keplerian fit for the combined Keck and Lick data sets of 
Marcy et al. (2001) has a much larger (22 ms -1 ) rms velocity scatter, and thus provides 
a less attractive starting point. This has motivated us to find a better analytic model for 
interacting systems which can bridge the gap between the dual-Keplerian approximation 
and the full three-body motion. 

We begin with a dual-Keplerian model using Jacobi coordinates. In this model, the 
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inner planet moves on an orbit about the star and the outer planet moves on an orbit 
around the center of mass of the inner two bodies. We then make the assumption that 
the two planets are undergoing librations about the 2:1 mean motion resonance. Marcy et 
al. (2001) have found that systems locked in this resonance can be stable on timescales 
of at least 70 million years, while other systems tend to be unstable. In our improved 
analytic model, the semi-major axes a of the two planets undergo sinusoidal oscillations 
about the exact resonant value in antiphase to each other. The period P res, amplitude A, 
and initial phase of these oscillations are treated as model parameters, in addition to the 
mean semi-major axis a 2 of the outer planet. The average mean motions fii and n 2 of the 
two planets are related by the condition that the rate of change of the resonance critical 
argument is zero at exact resonance. Hence 

fi\ — 2n 2 "b 2(712 ^ 1 ) = 0 (^) 


where 7 Ti and 7 r 2 are the longitudes of periastron of each planet. The average semi-major 
axes of the planets are related to their average mean motions by 

'G(m 0 + mi ) 1 1/3 


CLi = 


a 2 = 


ni 


G(m 0 + mi + ro 2 )] 1/3 


Tin 


( 2 ) 


where mo, mi and m 2 are the masses of the star, the inner and outer planets respectively. 
At a time t, the semi-major axes are given by 


Oi = ai[l + A X COS Tlre^t- t 0 ff)] 

a 2 = a 2 [l - A 2 cos n res (t - toff)] ( 3 ) 


where n res = 27r/P res and t 0 ff determines the initial phase of the resonant oscillations. In 
addition, conservation of energy requires that 

\(m 0 + mi)miaij 
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The initial mean anomaly Mv(Q) of the outer planet is treated as a model parameter, 
and the initial mean anomaly of the inner planet Mi(0) is then given by the critical 
argument for the resonance 

<7 = Ml — 2M2 + 2(7Tj — ^ 2 ) (5) 

where a = 0 at t = t oR . At time t, the mean anomaly of body i is given by 



The integrals are straightforward to evaluate since for each planet a, and hence n, is an 
analytic function of time. 

In this system, the mutual planetary perturbations are sufficiently strong that the 
longitudes of periastron will precess rapidly. We model this by allowing each periastron 
longitude to vary linearly with time, where the rates of change of the two angles represent 
additional free parameters. In principle, these parameters, and the resonance libration 
period P res , can be used to test the accuracy of the analytic model by comparing the 
precession rates with those from a full N-body integration. The orbits of the planets are 
assumed to be coplanar in our model, but the inclination i of this plane to the line of sight 
is included as a parameter. 

We used the analytic model to generate synthetic radial velocities for the central 
star and compared these with the observations from the Keck and Lick telescopes given 
in Marcy et al. (2001). We initially generated a randomized population of sets of model 
parameters and then used a genetic algorithm to evolve promising sets towards an improved 
description of the system. At each generation, the genetic algorithm evaluates the degree of 
fit for each parameter set, and cross breeds the best members of the population to produce 
a new generation. 

Figure 4 shows a model fit generated by the genetic algorithm for the Keck data alone. 
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This fit has a rms scatter of 8.1 ms \ which is comparable to the best dual-Keplerian fit 

V, 

or the preliminary fit obtained by the Levenberg-Marquardt N-body technique. Figure 5 
shows a model fit for the combined Keck and Lick data set. The rms scatter in this case is 
11.9 ms -1 , which represents a substantial improvement on the best dual-Keplerian model 
in Table 1. The apocentric orbital elements for the fits generated by the analytic model are 
given in Table 3. We again stress that other solutions with low rms scatter are likely to 
exist, and these orbital parameters, while being suggestive, do not necessarily represent the 
true dynamics of the system. 


5. Discussion 

The most important benefit of self-consistent dynamical fitting techniques for 
multi-planet systems is the ability to break the sini degeneracy and determine the true 
masses of the extrasolar planets. The true masses can occasionally be found in cases where 
the planet transits the parent star (e.g. Charbonneau et al 2000, Henry et al 2000), but 
such cases are unusual, and will be confined largely to planets with short periods. The 
foregoing techniques can in principle be applied to any system containing more than one 
planet, given a sufficient baseline of observation. Fischer et al 2000 have shown that roughly 
half of the planetary systems uncovered in the Lick Radial survey show evidence of a 
second companion. Thus we expect that numerous additional multi-planet systems will be 
forthcoming. Systems having massive short-period planets are especially amenable to this 
technique. As a specific example, an N-body integration of the Upsilon Andromedae system 
indicates that the planetary interactions are already producing observable deviations from 
the multiple Keplerian fit. 

Wolszczan (1994) has used a roughly similar analysis to the one described here to 
determine the true masses and inclinations of the planets orbiting pulsar PSR B1257+12. 
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However, in that case, a superposition of Keplerian fits provides a very good approximation 
to the observed reflex velocity of the pulsar. This stands in marked contrast to cases such 
as GJ 876, where the planetary interactions are an integral component of the overall motion 
of the star, and an analysis based on small perturbations to Keplerian motions may not 
necessarily succeed. 

There are several avenues for immediate additional improvement of the dynamical 
description of the emerging multi-planet exosolar systems. For cases where the radial 
velocity data is inadequate to cleanly delineate the planetary masses and orbital parameters, 
numerical integrations can reveal dynamical instabilities which can eliminate large portions 
of parameter space from consideration. The semi-analytic model of §4 can also be 
substantially improved through inclusion of realistic expressions for the secular evolution of 
the orbital elements (see e.g. Dermott & Murray 1999). 

We envision a four-stage procedure for determining the dynamical characteristics of 
multiple-planet systems in general. (1) As outlined by Butler et al 1996, periodograms (e.g. 
Lomb, 1976; Scargle 1982) can first be used to determine the fundamental periods within 
the system. (2) The radial velocity data can then be fit with fixed Keplerian ellipses (e.g. 
Butler et al 1999). (3) These multiple-Keplerian fits can then be further refined by versions 
of the semi-analytic scheme described in §4., and then (4) given a final self-consistent polish 
using the Levenberg-Marquardt scheme driving full N-body integrations. 
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Figure Captions 

Fig. 1 — Synthetic radial velocity variations for GJ 876 assuming a superposition of 2 fixed 
Keplerian motions for planets with elements given in Table 1 (red line), and an N-body 
integration using the same elements (black line). 

Fig. 2. — Synthetic radial velocity variations produced by an N-body integration using the 
osculating elements from Table 1 (Keck fit) and a best fit epoch of JD 2450671.98. 

Fig. 3. — Synthetic radial velocity variations for adit derived using the Levenberg-Marquardt 
procedure described in §3, using the Keck radial velocities. 

Fig. 4. — Synthetic radial velocity variations for a fit using the semi-analytic model described 
in §4, using the Keck radial velocities (solid circles). 

Fig. 5. — Synthetic radial velocity variations for a fit using the semi-analytic model described 
in §4, using the combined Keck (solid circles) & Lick (open circles) radial velocities. 



Parameter 


Inner Outer Inner 


Outer 


(Keck) (Keck & Lick) 


Period (day) 

30.1 

K (ms" 1 ) 

81 

Eccentricity 

0.11 

u (deg) 

328 

Periastron Time (JD) 



61.0 

30.12 

61.02 

211 

81 

210 

0.29 

0.27 

0.10 

329 

330 

333 


2450031.4 

2450106.2 


Table 1: Best-fit dual-Keplerian elements for GJ876 (from Marcy et al. 2001) 
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Parameter 

Inner 
*\ 4 * 

Outer 

Period (day) 

30.13 

61.58 

K (ms -1 ) 

80.9 

203.6 

Eccentricity 

0.226 

0.025 

u (deg) 

156 

70 

Mean anomaly (deg) 

277 

31 

a (AU) 

0.1297 

0.2092 

Epoch (JD) 

2450602.0931 


Table 2: Osculating elements derived by Levenberg-Marquardt N-body integra- 
tion scheme. 


Parameter 


Inner Outer Inner Outer 


(Keck) (Keck & Lick) 


Mass (Mj) 

0.740 

6.73 

1.927 

5.81 

a (AU) 

0.1302 

0.2077 

0.1298 

0.2082 

Eccentricity 

0.429 

0.240 

0.229 

0.006 

u> (deg) 

207 

179 

204 

97 

Mean anomaly (deg) 

186 

301 

136 

354 

Epoch (JD) 

2449995.0 

2449990.0 

sini 

3.34 

3.14 


Table 3: Osculating apocentric elements derived by the analytic model plus genetic 
algorithm scheme. 
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